library(stargazer)
library(plyr)
library(stringr)

load('stacked.frame.RData')

fit.w.robust <- function(model, cluster.var, dta){	
	robust.se <- function(model, cluster){
		require(sandwich)
		require(lmtest)
		M <- length(unique(cluster))
		N <- length(cluster)
		K <- model$rank
		dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
		uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
		rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
		rcse.se <- coeftest(model, rcse.cov)
		return(list(rcse.cov, rcse.se))
	}
	
	out <- model$model
	clustervar <- mapply(paste,cluster.var,dta[!(1:dim(dta)[1] %in% model$na.action),cluster.var],sep="")
	vcov <- robust.se(model, clustervar)
	model$se <- vcov[[1]]
	model$coeftest <- vcov[[2]]
	return(model)
}

stacked.model <- lm( reporters ~ treat*I(year.log==-3) + treat*I(year.log==-2) + treat*I(year.log==-1) + treat*I(year.log==1) + treat*I(year.log==2) + treat*I(year.log==3) + factor(paper.stack.fe) + factor(year.stack.fe) + circulation + home.county.circulation.share + population + white.share + median.income + bachelors.degree,data=stacked.frame)
stacked.model <- fit.w.robust(stacked.model,'paper.stack.fe',stacked.frame)

interactive.coefficients <- stacked.model$coefficients[grep(names(stacked.model$coefficients),pattern="treat:I(year.log",fixed=TRUE)]
interactive.ses <- stacked.model$coeftest[grep(names(stacked.model$coeftest[,2]),pattern="treat:I(year.log",fixed=TRUE),2]
interactive.upper <- interactive.coefficients+1.96*interactive.ses
interactive.lower <- interactive.coefficients-1.96*interactive.ses

plot.coefficients <- c(interactive.coefficients[1:3],0,interactive.coefficients[4:6])
plot.upper <- c(interactive.upper[1:3],0,interactive.upper[4:6])
plot.lower <- c(interactive.lower[1:3],0,interactive.lower[4:6])

pdf(height=4,width=6,file='figure2.pdf')
par(mar=c(4,5,.5,.5))
plot(x=c(-3,-2,-1,0,1,2,3),y=plot.coefficients,xlab="Years to Investment Purchase",ylim=c(-15,10),pch=16,las=1,type='n',ylab="Difference in Newsroom Employment
(Investment - Other Owners)")
abline(h=0,col='gray85',lwd=3)
abline(v=0,lty=2)
points(x=c(-3,-2,-1,0,1,2,3),y=plot.coefficients,pch=16,cex=1.5)
lines(x=c(-3,-2,-1,0,1,2,3),y=plot.coefficients,lwd=3)
segments(x0=c(-3,-2,-1,1,2,3),y0=interactive.upper,y1=interactive.lower,lwd=3)
dev.off()

stacked.model.entertainment <- lm( entertainment.reporters ~ treat*I(year.log==-3) + treat*I(year.log==-2) + treat*I(year.log==-1) + treat*I(year.log==1) + treat*I(year.log==2) + treat*I(year.log==3) + factor(paper.stack.fe) + factor(year.stack.fe) + circulation + home.county.circulation.share + population + white.share + median.income + bachelors.degree,data=stacked.frame)
stacked.model.entertainment <- fit.w.robust(stacked.model.entertainment,'actual.name',stacked.frame)

interactive.coefficients.entertainment <- stacked.model.entertainment$coefficients[grep(names(stacked.model.entertainment$coefficients),pattern="treat:I(year.log",fixed=TRUE)]
interactive.ses.entertainment <- stacked.model.entertainment$coeftest[grep(names(stacked.model.entertainment$coeftest[,2]),pattern="treat:I(year.log",fixed=TRUE),2]
interactive.upper.entertainment <- interactive.coefficients.entertainment+1.96*interactive.ses.entertainment
interactive.lower.entertainment <- interactive.coefficients.entertainment-1.96*interactive.ses.entertainment

plot.coefficients.entertainment <- c(interactive.coefficients.entertainment[1:3],0,interactive.coefficients.entertainment[4:6])
plot.upper.entertainment <- c(interactive.upper.entertainment[1:3],0,interactive.upper.entertainment[4:6])
plot.lower.entertainment <- c(interactive.lower.entertainment[1:3],0,interactive.lower.entertainment[4:6])

stacked.model.other <- lm( I(reporters-entertainment.reporters) ~ treat*I(year.log==-3) + treat*I(year.log==-2) + treat*I(year.log==-1) + treat*I(year.log==1) + treat*I(year.log==2) + treat*I(year.log==3) + factor(paper.stack.fe) + factor(year.stack.fe) + circulation + home.county.circulation.share + population + white.share + median.income + bachelors.degree,data=stacked.frame)
stacked.model.other <- fit.w.robust(stacked.model.other,'actual.name',stacked.frame)

interactive.coefficients.other <- stacked.model.other$coefficients[grep(names(stacked.model.other$coefficients),pattern="treat:I(year.log",fixed=TRUE)]
interactive.ses.other <- stacked.model.other$coeftest[grep(names(stacked.model.other$coeftest[,2]),pattern="treat:I(year.log",fixed=TRUE),2]
interactive.upper.other <- interactive.coefficients.other+1.96*interactive.ses.other
interactive.lower.other <- interactive.coefficients.other-1.96*interactive.ses.other

plot.coefficients.other <- c(interactive.coefficients.other[1:3],0,interactive.coefficients.other[4:6])
plot.upper.other <- c(interactive.upper.other[1:3],0,interactive.upper.other[4:6])
plot.lower.other <- c(interactive.lower.other[1:3],0,interactive.lower.other[4:6])

x.seq <- c(-3-.075,-2-.075,-1-.075,0,1-.075,2-.075,3-.075)
x.seq.2 <- c(-3+.075,-2+.075,-1+.075,0,1+.075,2+.075,3+.075)

pdf(height=4,width=6,file='figure3.pdf')
par(mar=c(4,5,.5,.5))
plot(x=c(-3,-2,-1,0,1,2,3),y=plot.coefficients.other,xlab="Years to Investment Purchase",ylim=c(-15,10),pch=16,las=1,type='n',ylab="Difference in Newsroom Employment
(Investment - Other Owners)")
abline(h=0,col='gray85',lwd=3)
abline(v=0,lty=2)

points(x=x.seq,y=plot.coefficients.entertainment,pch=16,cex=1.5,col='gray50')
lines(x=x.seq,y=plot.coefficients.entertainment,lwd=3,col='gray50')
segments(x0=x.seq,y0=plot.upper.entertainment,y1=plot.lower.entertainment,lwd=3,col='gray50')

points(x=x.seq.2,y=plot.coefficients.other,pch=16,cex=1.5)
lines(x=x.seq.2,y=plot.coefficients.other,lwd=3)
segments(x0=x.seq.2,x1=x.seq.2,y0=plot.upper.other,y1=plot.lower.other,lwd=3)

legend("topright",c('General/Political','Entertainment'),pch=c(16,16),col=c('black','gray50'),cex=1,pt.cex=1.5)

dev.off()
